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Abstract 

The effectiveness of the recently developed Fixed-Node Quantum Monte 
Carlo method for lattice fermions, developed by van Leeuwen and co-workers, 
is tested by applying it to the Id Kondo lattice, an example of a one- 
dimensional model with a sign problem. The principles of this method and 
its implementation for the Kondo Lattice Model are discussed in detail. We 
compare the fixed-node upper bound for the ground state energy at half filling 
with exact-diagonalization results from the literature, and determine several 
spin correlation functions. Our 'best estimates' for the ground state correla- 
tion functions do not depend sensitively on the input trial wave function of 
the fixed-node projection, and are reasonably close to the exact values. We 
also calculate the spin gap of the model with the Fixed-Node Monte Carlo 
method. For this it is necessary to use a many-Slater-determinant trial state. 
The lowest-energy spin excitation is a running spin soliton with wave number 
7T, in agreement with earlier calculations. 

dedicated to Hans van Leeuwen at the occasion of his 65 th birthday 
PACS numbers: 71.27.+a,71.10.+x, 75.10.Jm,71.20.Ad 
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I. INTRODUCTION 



As is well known, quantum Monte Carlo simulations are plagued by the so-called sign 
problemSS. The sign problem refers to the fact that when physical properties are sampled 
in configurations space, one collects large positive and negative contributions due to the fact 
that a fermion wavefunction is of different sign in different regions of configuration space. 
These contributions of opposite sign tend to cancel, giving results that may be exponentially 
smaller than the separate positive and negative contributions. Though the sign problem can 
be circumvented in special cases, e.g., for the Hubbard model at half filling!, no general 
solution has emerged yet from the various approaches that have been explored to cure 
iti~0. 

In 1990, when B. J. Alder was Lorentz Professor in Leiden, Hans van Leeuwen became 
acquinted with the Fixed Node Monte Carlo (FNMC) method of Ceperley and AlderEMi, 
which avoids the sign problem in the context of continuum Green's function Monte Carlo. 
This stimulated him to explore the possibility of formulating a lattice version of FNMC, first 
with a postdoc, Ani, and later in collaborations with the present authors00. The formula- 
tion of the approach which was developed later0 was shown to be variational^, i.e. to give 
an upper bound to the exact ground state, and is the subject of this paper, which we dedicate 
to Hans van Leeuwen. We test this FNMC method for lattice fermionsSfllll on a simple one- 
dimensional (Id) model for which various results are available, the Id Kondo lattice model 
(KLM) at half filling. This FNMC method involves an approximation that removes the sign 
problem in the context of Green's function Monte Carlo. Different Monte Carlo techniques 
that have been applied to the Id KLM include the world-line algorithm^, a finite temper- 
ature grand-canonical method involving a Hubbard-Stratonovich transformationS and the 
ground state method developed by Sorella et aZ.0'El. All three Monte Carlo methods suffer 
from the sign problem, even in Id. 

The KLM is one of the basic models for correlated fermions. It can be obtained as the 
strong-coupling limit of the periodic Anderson model, which aims at capturing the essen- 
tial physics of heavy-fermion materialsiiB. In the limit of strong on-site repulsion among 
the /-electrons, a picture emerges of localized /-electrons interacting with a conduction 
band. In recent years, the model has been studied by a variety of methods, including varia- 
tional approximations, exact diagonalizations and the density matrix renormaiization group 
met ho di §S@H This, together with the fact that quantum Monte Carlo simulations of 
this model do have a sign problem, makes the Id KLM a suitable testing ground for our 
lattice FNMC method. 

The lattice version of FNMC gives, like the continuum version0Hll0 which inspired it, 
upper bounds for the ground state energy01l. It improves upon a trial wave function for 
a given Hamiltonian by employing a Green's function projection method with a modified 
Hamiltonian in which all terms which would lead to unwanted sign changes in the sampling, 
are treated in a special way. The sign structure of the resulting approximate wave function, 
which is the ground state wave function of the modified Hamiltonian, is the same as that of 
the original trial wave. One obvious immediate question of interest is how close the FNMC 
energy estimate is to the exact ground state energy of the original Hamiltonian. We will 
study this question by applying the fixed-node projection to a trial wave function with a free 
parameter. As we shall see, for the KLM, the ground state energy obtained in the FNMC 
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at half filling is quite independent of the precise trial wave, and quite close to the values 
obtained in exact diagonalisation. We also compare the FNMC results for some correlation 
functions (which do not obey bounds), with exact result for chains of six sites, with 

coupling constant J equal to 0.2 and 1.0. Here too we find that our best FNMC simulation 
estimates are rather independent of the starting trial wave function, and reasonably close to 
the exact values. We finally also show how the spin gap of the Id KLM can be determined 
with our FNMC, although this is computationally much more demanding, since a trial wave 
consisting of a sum of slater determinants must be used. Good agreement is found with 
earlier results§H in this case. 

Before presenting our results, we first briefly discuss the Id KLM and the reason why 
sampling it with unrestricted random walks leads to the sign problem. Then, in section |T|, 
the principles of the lattice FNMC are summarized, followed by details of the implementation 
for the KLM. Section [IV] gives the comparison with exact results for small lattices. In section 
[V], our results for the running spin triplet excitation are discussed. 



II. THE KONDO LATTICE MODEL 

The Kondo lattice Hamiltonian is given by 

Kklm = -t ( c L c jv + c )a c i<j) ~ l^^2n ic + J^2S if ■ S ic ■ (1) 

(ij) i i 

The two kinds of electrons, denoted by c for the conduction band and / for localized levels, 
have a spin-spin interaction. For the /'s, the constraint is that there has to be precisely one 
/-electron on every site. In (JJ) and below, a summation convection is used for repeated spin 
indices. 

We seek to write the Hamiltonian in a form that is convenient for GFMC calculations. 
If one would treat the /'s as spins, which are not antisymmetrized but form a dynamical 
background for the conduction electrons, the total number of up-spin (and of down-spin) 
fermions would not be conserved by the Hamiltonian. It is not convenient to use this 
representation in a GFMC calculation whose starting trial wave is a fully fermionic mean- 
field type wavefunction. If we use the constraint of one /-electron per site and the identities 

I -> 1 4. 

Sif — —fivTaa'fia' j Sic = ~Z C ia^cr(T' C ia' i (2) 

where the components of T au i denote the three Pauli matrices, the Hamiltonian can be 
written in the form 

WKLM = -tJ2 {4cr C j<? + c]*Ci*) ~^Y1 ( C lfi°) {fia'Cia') ~ n' Yl U ^ , (3) 

(ij) i i 

with fi'—fi— J/4. This is now fully in fermionic language. Therefore, Slater determinants of 
fixed dimension can be used. 
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FIG. 1. A sequence of processes that lets two conduction electrons of equal spin pass each other in Id. 
Large arrows denote /-electrons, small ones denote c-electrons. The successive states in the sequence are 
separated by vertical dashed lines. 

In many fermionic lattice models, e.g. the Hubbard model, the fermion statistics is not 
really important in Monte Carlo simulations in Id. The reason is that, fermions of the same 
spin cannot pass each other in Id and so their ordering is fixed. Since the exchanges that 
give rise to sign changes as a result of the fermion antisymmetry of the wave function, are 
suppressed, there is no sign problem. In higher dimensions, this is not the case. 

For the Id KLM, there is no fixed ordering of conduction electrons of equal spin. The 
presence of the spin flip term proportional to J in Hklm makes it possible for the ordering 
to change. This is illustrated in Fig. 1. A down c-electron at site 1 first changes in an 
up-electron, due to the simultaneous flip of a c-f pair. Then, the c-electrons at site 1 and 2 
have opposite spins, so after a hop due to the kinetic term, they can occupy the same site. 
After an additional hop of the other c-electron, the two c-electrons then effectively have been 
interchanged. 

Such an interchange can also happen in a Monte Carlo simulation, and one needs to 
take into account that two configurations that differ by the interchange of two numbered 
electrons must have opposite signs in the wave function. This is the reason why even the Id 
KLM exhibits a sign DroblemHEl. In the case of 6 sites at half-filling with J = 1.0, which 
is studied by Otsuka", it appears that the sign problem is not very severe, but at certain 
filling fractions the sign problem is known to make simulations prohibitively difficult!!. 

The Id KLM has been studied in different regimes. If the number of /-electrons is eoual 
to the number of sites, and the carrier concentration is low, there is a ferromagnetic state@H. 
In weak coupling, at larger density but below half-filling, one obtains a paramagnetic stated, 
and at half-filling, the system shows insulating spin-liquid behaviorSH. Recently, the ground 
state was proven to bea spin singlet andproven to be uniqueii. For slightly less than one 
/-electron per site, impurity bands ariseEa 

In this paper, we limit ourselves to the case of half-filling, one c-electron per site. Finite- 
size scaling resultsll as well as recent density matrix renormalization group calculations!! 
both show that there is a gap for spin and charge excitations for all J > and thus confirm 



4 



the insulating spin-liquid character of the ground state at half filling. 



III. THE FNMC METHOD FOR THE KLM 



A. Principles of the FNMC method 



Since the general principles of the FNMC for lattice fermions have been laid out 
beforeM-acI, we only summarize the most essential aspects of the method here. 

In a Green's function Monte Carlo method, one projects out the jrround state of a system 
with Hamiltonian Ti from an initial trial state \iPt)- As beforeSSi, we use a projection 
operator T which acts as follows: 

\r)=T n \^) = \\-T{H-w)\ n \^). (4) 

The implementation is in configuration space, and has a stochastic character. A specific 
configuration in configuration space, which determines the locations, spins, etc. of all the 
labeled electrons, is denoted by R, and we write ^t(-R) = (Rh^r), etc. When r is small 
enough and w adjusted properly during the sampling processi'O, the operator T n projects 
onto the ground state as n — > oo. To obtain better statistics, we introduce importance 
sampling: we let the Green's function 

G(R, R') = ip T (R)F(R, R')^ X {R') (5) 

determine the transition probabilities of a random walker from R' to R; for simplicity, we 
take the trial wave function to be real. The projection (§) then becomes 

r{R) = E MR)G(R, R n )G(R n , iVi) ■ • ■ G(R 3 , R 2 )G(R 2 , RM 2 T {R X ) (6) 

In the random walk interpretation underlying the Monte Carlo process, the initial distribu- 
tion of the random walkers is given by ip^(R), and is sampled stochastically by splitting 

G as 



G(R, R!) = P(R, R')m(R') , (7) 

with m(R') =J2rG(R, R') and hence J2rP(R, R') — 1- This notation anticipates that we 
wish to view P as a transition probability, the probability for a particle to make a transition 
from configuration R' to R, so that a path Ri, R 2 , ■ ■ ■ , R n in configuration space is generated 
by sampling the transitions according to P. The weight factors m which are accumulated 
along a path are sampled by viewing them as a multiplicity factor of each walker^ED. After 
a suitable number of steps, these multiplicity factors are sampled by a branching process: 
at these events, a walker with a multiplicity factor m can be either killed, stay alive, or split 
into more walkers in such a way that, on average, there are <m> new ones after the event. 
After each branching event, the factors m of all the walkers are reset to 1. 

If all transition probabilities G would be positive, the above process could be imple- 
mented straightforwardly as in simulations of boson lattice models^. For fermions, the 
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sign problem arises in the present formulation through the fact that G(R,R'), and hence 
P(R,R'), can be negative. In particular, for the KLM (|3|), as for the Hubbard mo defi all 
off-diagonal terms (R\H KLM \R') of the Hamiltonian TLklm are negative, and so negative 
signs arise in making transitions between configurations at which the trial wave function has 
opposite signs, x/jt(R)/iPt(R') < 0. Taking those transition probabilities P(R,R') propor- 
tional to \G(R, R')\ and carrying the sign with each walker would give positive and negative 
multiplicities, eventually, and therefore, the implementation of Eq.(^) would lead to large 
positive and negative contributions to all measured quantities. The resulting cancellations 
are the essence of the sign problem. 

In the lattice FNMC the sign problem is avoided by introducing a slightly different 
effective Hamiltonian H. e ff in which all steps that lead to a sign change of the walker are 
left out, and replaced by a potential term0'3. The steps that are left out satisfy 

{R\H\R'}ijj T (R)ip T (R') > 0. (8) 

This is implemented by defining an effective Hamiltonian as follows: The off-diagonal terms 
are 

(R\H cS \R'} = (R\H\R') (if (R\H\R')iIj t (R)^ t (R') < 0) , 

= (otherwise), (9) 

and the diagonal terms are 

(R\H cS \R) = (R\H\R) + (R\V S{ \R), (10) 
where the 'sign flip' potential that replaces the hops that satisfy Eq.(|S|) is given by 

(R\ Vsl \R) = jr(R\H\R')^^. (11) 



R> 



MR) ' 



In this expression, the sum runs over all configurations R' connected by a non-zero matrix 
element (R'\H\R) to the configurations R, for which (|8|) holds. The ground state energy of 
Ti-eff-, which can be sampled without sign problem, give sS an upper bound to the ground 
state energy of the true Hamiltonian Tt. Expectation values of physical quantities are then 
obtained in the standard 



B. The variational mean-field type trial state for the KLM 

As we saw above, a prerequisite for a Green's function Monte Carlo calculation is a trial 
wave function. For the KLM, we use what amounts to a Gutzwiller-projected mean-field 
type wavefunction as a trial state. This wavefunction is essentially obtained as follows. Since 
earlier work indicates that this gives the lowest energy results, we use the Kondo decoupling 
scheme in Eq. (§), i.e., TLklm is approximated by 

H v = -tJ2 {clc ja + Cj W) - £ [(cifj) V t + (flcj) V*] + -Y, \V t \ 2 . (12) 



(ij) 
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The Hamiltonian Tiy is bilinear in the fermion operators and hence can be diagonalized. 
We denote the ground state of this Hamiltonian by |Vv)- The trial state \ipr) we use for our 
calculations is then 

|V, T ) = P G |W>, (13) 

where Pq is the Gutzwiller projection operator which projects onto states in which each site 
is occupied by one /-electronc3EJ. 

For the homogeneous ground state, all V^'s should be taken equal, Vi — V. Following 
Otsuka§, we will use this V as a variational parameter to construct a family of trial states for 
our ground state calculations. The explicit form of the wavefunction can be easily obtained, 
and is given by Eq. (5) of OtsukaS (his parameter V is the same as ours) . This wavefunction 
takes the form of a hybridized band state, but after Gutzwiller projection, it can also be 
written in the form of a overlapping Kondo cloud stated. 

We will actually find it more convenient to present our ground state results as a function 

of 

b = ^ v \flc ia \^ v ). (14) 

Note that the average is computed before Gutzwiller projection - in the Gutzwiller projected 
state the average is obviously zero. In the mean-field approximation, the self-consistency 
condition for the homogeneous ground state reads Jb/2 = V; this relation can easily be 
worked out in the thermodynamics limit, but as stated before, we will not use this. 

We also use the Gutzwiller-projected mean field solution as our trial state for the lowest 
energy triplet state in section V. The mean-field solution in this case is inhomogeneous0; 
hence, in this case the parameters Vi in the selfconsistency conditions J/2('0vJ// a .Cj (r |'0v) = Vi 
do depend on the spatial index i. We refer to the paper by Wang et a/£il for a detailed 
discussion of the structure of this mean field solution. 

The single-particle levels or our trial wave are represented by an index for the energy 
the level, an index for the site, and an index which indicates whether an electron has cor / 
character. This way of representing the trial wave function is suitable for the order in which 
the operators in the interaction term appear in Eq.(0), and for the decoupling we have chosen 
to generate the trial state. The operators between parentheses in the spin interaction term 
in TLklm represent intermediate steps in a Monte Carlo diffusion process. These correspond 
to changes within one spin sector. It is, therefore, natural to have numbered electrons of a 
certain spin and to allow changes from c to / and vice versa, rather than to have numbered 
electrons with the c or / character fixed and letting the spin change. Both representations 
are equivalent, but our choice allows to work with Slater- matrices of fixed size. 

In the Monte Carlo calculation, the trial wave determines the distribution of random 
walkers; each walker represents a configuration, i.e., specifies the positions of each electron, 
its spin, and whether it is c or /. The weight of a certain configuration in the initial 
ensemble can be calculated from the trial state, by taking the product of the determinants 
corresponding to the spin-up and spin-down single-particle states. The ensemble is chosen by 
generating configurations at random, and then comparing the weight squared with a random 
number, in order to decide whether that configuration should be accepted as a member of 
the ensemble or not. By imposing the constraint of one /-electron per site, the ensemble is 
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automatically Gutzwiller-projected. Once the initial ensemble is Gutzwiller projected, the 
ensemble remains so during the projection process, since all moves allowed by TLklm keep 
the /-levels singly occupied. 



C. Implementation of the FNMC for the KLM 

For the KLM Hamiltonian Eq.(|3|), all off-diagonal terms (R\TL\R) are negative, for an 
antiferromagnetic spin-interaction J>0. All steps are therefore according to (J8|) subdivided 
into allowed steps for which ip T (R)ip T (R') > and forbidden steps for which x(j T (R)ip T (R') < 0; 
the latter steps contribute to the sign flip potentialil (|TTD. If we use TLklm in the projector 
(U), three things can happen in one time step of the FNMC. First of all, R' can go to a 
configuration with one c-electron hopped to a neighboring position due to the kinetic term 
in TLklm- The second possibility is a simultaneous spin flip: the spin interaction term 
proportional to J allows a configuration which has, on a certain site, a pair c-up, /-down, 
to change into c-down, /-up (or vice versa). The third possibility is that nothing happens 
in a time step: R' = R' . The relative probabilities are given by Eq. (|5[). For a given 
walker, which corresponds to a given configuration, a list is therefore made of all possible 
allowed steps and their probabilies. When a forbidden step is encountered in making this list, 
the corresponding contribution to the sign flip potential ( pTT| ) is calculated. Since for every 
configuration R at most one electron changes its state (the site- or c/f label) per spin sector, 
the ratio ipr(R) /iPt(R') which determines the probabilities and V s {, can be calculated in a 
number of operations that is linear in the size of the system, if one already has the transposed 
inverse of the Slater matrices availableil. 

Once an ensemble of random walkers with weight determined by the trial wave function 
has been prepared, as described in the previous subsection, the Monte Carlo projection 
is done according to For a given walker, all possible moves are considered, and for 
each move, the ratio's ^t(R) j 'iPt(R') are calculated. This operation corresponds to a dot- 
produdU, so the time needed to compute it is linear in the system size. If ijjT(R)/ijjT(R') 
is positive, the step to R is allowed. The probability factor G(R, R') = P(R, R')m(R') 
of allowed Monte Carlo moves are stored in an ordered table in which each element is 
the sum of the previous element and the probability factor r\ipr{R) /Vt(-R')I f° r a hop or 
Jt\i/jt(R)/i/jt(R')\/2 for a spin flip. The last element is the sum of the one but last element 
and the probability factor for staying, G(R,R') = 1 — r(U pot + V s { — w), where U pot + V s f 
is the total potential energy of TL e ff. In this way, the value of the last element equals 
J2rG(R, R') = m(R'); the random decision to select a move or to stay is then made by 
deciding between which elements of the ordered table the product of a random number 
between and 1 and m(R') falls, using the Numerical Recipes routine located. 

The first stage of the diffusion implementation of the projection is a thermalization. 
During this stage, the parameter w in (HJ) T = 1 — r{TL — w) is adjusted in such a way that 
the ensemble of walkers stabilizes under the branching process by which the multiplicities 
of the walkers are updated; w approaches the measured ground state energy in this process. 
After the thermalization, the usual quantity being measured in a Green's function Monte 
Carlo is the mixed estimator (i/)x\0\ip n ) for the local value of an operator O, 
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O local (R) = (R\0\xh)/MR) 



J2(R\0\R')(R'\^) 

R' 



/MR)- (is) 



The mixed-estimator is directly measured in the FNMC program, but a better estimate for 
an expectation value is obtained by assuming that the trial state is close to the ground state 
and neglecting quadratic terms in the difference^: 

<^ n | W") « 2 • (MOm - (ih\0\,fa). (16) 

We will refer to the right hand side as the best estimate. 

Operators which are diagonal, like the S z spin correlation function are most efficiently 
calculated, since for off-diagonal terms computation of the ratio's iPt(R')/4>t(R) takes sub- 
stantial computer time. In our FNMC, r needs to be small enough that T n projects onto 
the ground state, and large enough that the convergence is sufficiently rapid§. We typically 
work with an ensemble of on average N ens = 1000 walkers. In one interval, all walkers are 
propagated during N time time steps before branching. iV time is chosen such that the multi- 
plicity factors m remain less than 2. After a thermalization of iVtherm intervals, statistics 
is accumulated in iVbi oc k blocks of A^ ntv intervals each. In principle, the blocks are treated 
as independent measurements, and occasionally we check whether these are sufficiently in- 
dependent indeed. If necessary, we increase iVi n t v to make them more independent; an 
example of this will be discussed in section V. The values of all these parameters used in 
the simulations will be listed in the figure captions. 



IV. RESULTS FOR J = 0.2 AND J = 1.0 

As a first test of the FNMC on the KLM we compare with exact diagonalization results 
by Yamamoto and Uedai. The coupling constant is J = 0.2, and the system consists of six 
sites with periodic boundary conditions. The trial wave functions we use are as described 
in section IV. B. 



In Fig. 2 we plot the FNMC energy as a function of 6, defined in (14]V. Note that the 
energy estimates are above the exact ground state energy, as they should beE-i Furthermore, 
we see that the minimum is quite flat in the range 0.15<5<0.7, and very close to the exact 
value (note the vertical scale!). Thus, our estimate for the ground state energy one is quite 
independent of the trial wave function — apparently, therefore, while the variational energies 
do depend stronly on b, the projected energies are not very much affected by the fixed-node 
constraints over some range of values of b. Finally, also note that the statistical fluctuations 
are smaller close to the optimal value of b, 6^0.25, in agreement with the general trend that 
fluctuations are smaller the better the ground state is approximated, and that statistical 
fluctuations are reduced if there is a gap in the excitation spectrumcll (the Id KLM does 
have a gap). 
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FIG. 2. Energy of a six-site KLM with J = 0.2; the parameters used in the Monte Carlo calculation 
are r = 0.01, w stai t = -8.3, N timc = 30, N thcim = 20, N intv = 5, iV b iock = 400, iV cns = 1000. The dashed 
horizontal line denotes the exact result of Yamamoto and UedalEI. 
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FIG. 3. Two examples of correlation functions in the J = 0.2 KLM. The parameters used in the FNMC 
program are r = 0.03, w stalt = -8.3, N timc = 30, N thcim = 20, N intv = 5, N hlock = 200, = 1000. Only 
the mixed estimators are shown. The short dashed line indicates the exact value from Ref.EZI, and the long 
dashed lines indicate the precision to which this exact result was given. 
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An example of the mixed FNMC estimates of two spin correlation function, {iPt\S 1 c - S l c \tl) n ) 
and (V'tI'S'/ • Sl\ip n ) for J = 0.2, are shown in Fig. 3. We see that near the optimal value of 
b, the mixed estimate is close to the exact value. At the same time, this figure illustrates 
that, not unexpectedly, the correlation functions are more sensitively dependent on the trial 
wave function than the projected energy shown in Fig. 2: at b~0A, the mixed estimate of 
the nearest neigbor c-f spin correlation function is almost a factor of 2 off, while the energy 
at this value is still quite close to the proper value. As we shall illustrate in detail below for 
,1 — 1, the estimates improve when we consider the best estimates instead. 

For J =1.0 we follow the same procedure as for J = 0.2 and compare with exact results 
obtained by Otsukalij. The upper line in Fig. 4a gives the energy measured in the starting 
ensemble, the lower line is the Fixed-Node value, i.e., after projection. Like in the case 
of J = 0.2, the latter curve is very flat, while the starting values depend strongly on the 
input wave function. Fig. 4b shows the same data on an expanded scale, on which one 
can see that the flat part of Fig. 4a really has a minimum. The exact diagonalization 
result" E = —8.561616 is also indicated in the picture. Clearly, also in this case the FNMC 
projection method is able to come quite close to the exact energy even if we start from a 
trial state that has a bad energy, and the statistical fluctuations are again smallest close to 
the minimum in FNMC energy. 
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FIG. 4. Energy of a six-site KLM with J = 1.0; the parameters used in the FNMC calculation are 
r = 0.003, w start = -9.9, Mime - 20, JV thorm = 20, N intv = 1, AWk = 10000, N CDS = 1000. In the right 
plot, the same results are plotted on an expanded scale. 

In Figs. 5 and 6, we present the results for on-site correlation functions and correlation 
functions involving different sites, respectively. Three values are plotted: the variational 
value (using the Gutzwiller projected state |Vv}), the mixed estimator and the best estimate 
given by fliTf). For all correlation functions considered, the latter curve is relatively flat 
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throughout the whole range of b values where the projected energy shown in Fig. 4 is close 
to the exact value. Thus, although the estimated correlations are slightly off, these results 
do show that, at least for the KLM, correlation functions are not strongly dependent on the 
trial state, and hence can be estimated relatively well with our FNMC. 
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FIG. 5. On-site correlations in the J = 1.0 KLM. The parameters used in the FNMC calculation are 
r = 0.003, = -9.9, Mime = 20, iV therm = 20, N intv = 1, 7V block = 200, N cns = 1000. 
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FIG. 6. Correlations on different sites in the J = 1.0 KLM. The parameters used in the FNMC 
calculation are t = 0.003, w start = -9.9, iV timc = 20, AT therm = 20, N lntv = 1, AWk = 200, iVc„s = 1000. 

V. FNMC CALCULATION FOR THE SPIN SOLITON 

The lowest-energy excitation bove the S = ground state of the half-filled KLM has total 
spin S = l§fi In a mean-field calculation, one is able to obtain a self-consistent solution 
with the spin excitation localized on a few sites0. Wang et a/.0 proceed by performing a 
Gutzwiller-projected mean-field calculation and by writing 

= H ex P(^ x c)|VO, (17) 

with \ip Xc ) = 'Pg\'4 ) ™/) the Gutzwiller-projected local triplet state with the center of the 
soliton located at x c . The minimum of this dispersion is at wave number q = it. We follow 
the general strategy of investigating the robustness of mean-field results by using this wave 
function as trial wave function in a FNMC calculation. To obtain the spin-gap in the FNMC, 
we perform calculations both in the 5 = and in the 5=1 sector. GFMC does not always 
project on the ground state, only on the lowest state that has a component along the trial- 
state. Here, we use this to our advantage: the total spin is conserved by the Hamiltonian, 
and, therefore, if one starts in the S = 1 sector, one remains in the S = 1 sector. Comparing 
the lowest energies in both sectors gives the gap. 

Eq. ( |17D as it stands, seems to indicate that not only different signs occur, but also 
different complex phases. Because of reflection symmetry, however, one can combine q and 
— q and write 

l^ff) = Hcos(gx c )|^ Xc ), (18) 
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which is a real problem againu. 

We perform FNMC calculations for a system of 20 sites with periodic boundary condi- 
tions. This however takes computer time: for each possible step, 20 ratios of determinants 
need to be calculated, not just one, since 



The factors ip Xc (R) / ip Xc (R') can be obtained as simple dotproducts againt^l, in most cases. 
The exception is when i/; Xc (R') =0: in those cases ip Xc {R) needs to be calculated from scratch, 
for which the computer time increases as the cube of the size of the system. 

Note that this difficulty never occurs in a calculation with only one Slater determinant 
as trial wave function: with importance sampling, the probability to be in a configuration is 
proportional to iPt{R'), so if this is zero, a random walker never visits such a configuration. 
Therefore, we never need to compute ratio's iPt{R) / Vt(-R') in which the old configuration 
R' has ip Xc (R') — 0. In the present case, only J2 Xc cos{qx c ){R'\ip Xc ) determines the probability 
to visit a configuration R', and a single ip Xc (R') may be small or zero. 

While smallness may be a practical difficulty in terms of numerical accuracy, the main 
problem is that a Slater-matrix can be really singular. In the case of a S — 1 soliton trial- 
state with S z — 1, this turns out to happen for (R'\ip Xc ) if the /-electron at site Xr has spin 
down, in configuration R', as is illustrated by the results of Fig. 2 of Wang et a/ril. 

Considering possible moves, ip Xc (R) needs to be calculated for all x c and for all R. Com- 
putation of many values ipx c {R) from scratch would be very time-consuming. However, the 
only possibility to make a singular matrix non-singular, is to flip the spin of the /-electron 
(and of the c-electron) on site x c . Such a flip can make one singular matrix non-singular, 
and only for the corresponding new configurations R do we need to calculate one Slater 
determinant. 

In practice, we keep track of which of the 20 matrices are singular, for a certain random 
walker. For all hops, all non-zero new values ip Xc (R) can be calculated as dotproducts, so 
that in calculating the ratio fll~9|) each flip in which one goes from /-down to /-up can make 
one Slater matrix non-singular, and for this one the determinant has to be calculated from 
scratch. 

Once a step has been chosen, all matrices have to be updated. If it is a flip from /-up 
to /-down, one determinant becomes singular. If it is a flip with f-down to f-up, one matrix 
becomes non-singular, and its transposed inverse needs to be calculated, for facilitating the 
calculation of transition probabilities of subsequent steps. Except for this more complicated 
calculation of iPt(R) / Vt(-R') > the FNMC program is the same as described before. 
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FIG. 7. Illustration of how a run proceeds: The average Monte Carlo energy calculated during one 
'block' in the FNMC calculation for a k = tt soliton in a KLM on a Id lattice of 20 sites with J = 4.0, 
calculated with parameters t = 0.01, w start = -90, N therm = 3, N time = 1, N lntv = 3, N Uoc k = 400, 
N ens — 200. The inset shows the error bar of the last 300 blocks, and calculated by grouping N group 
measurements together. For N group > 3, the error bars do not increase, indicating that the energies of 
different groups are statistically independent. 

Fig. 7 illustrates how the FNMC projection proceeds as the number of iterations in- 
creases. First, one observes projection on the ground state: the energy measured over a 
number of Monte Carlo time steps, a 'block', drops. Then, there are fluctuations around a 
mean value. Accumulating statistics leads to a reduction of the error bars. One observes 
that the measurements in consecutive blocks are not independent. The correct error bar is 
obtained by grouping N group measurements together and thus dividing the iV tota i blocks in 
-^totai /-Wgroup measurements. After this regrouping, one has fewer values, but they are more 
independent. The error bars this gives are plotted in the inset of Fig. 7: the value of the 
plateau is the error bar we report in the energy dispersion curve (all error bars in energies 
are obtained this way). 
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FIG. 8. Dispersion of a running soliton in a J = 4.0 KLM of 20 sites. The parameters used in the 
Monte Carlo calculations are r = 0.005, u> s tart = —90, iVtime = 1, A^hcrm = 3, JV intv = 3, iVbi oc k = 200, 



The FNMC dispersion of the spin soliton, for J = 4.0, on 20 sites, is presented in Fig. 
8. Since the S = value is Egz$ = —63.423(5), and the minimum of the dispersion is 
at E^ MC = -60.47(3), the gap we obtain in FNMC approximation is A™ MC = 2.95(3), 
which is in agreement with the results shown in Fig. 3 of Wanget al. using a Gutzwiller- 
projected mean-field approximation, and those of Yu and WhiteH using the density matrix 
renormalization group method. So, while the energies in both the S — and the S = 1 sector 
drop relative to the one estimated with the Gutzwiller-projected mean-field wave function, 
the difference between the S = and the S = 1 ground state energies is essentially the same 
as what is known to be the correct valueS 



VI. CONCLUSIONS 



We have applied the lattice FNMC method to the Id KLM. We observe that, for small 
lattices, quite accurate ground state energy estimates are obtained, even if the starting trial 
wave function gives quite a bad approximation to the energy. In the FNMC, the energy 
estimate is always above the exact ground state energy, but in the cases studied here, they 
are only slightly above the exact values for a large range of values in the variational trial 
wave. Over the same range over which a good ground state energy estimate is obtained, 
reasonably accurate values for correlation functions in FNMC approximation are obtained, 
which again are rather independent of the starting trial-state of the KLM. To be able to 
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calculate a dispersion of an excitation, one needs a conserved quantity, such that the lowest 
energies in different sectors can be compared — e.g., the isospin gapEa of the KLM can not 
be obtained within our FNMC, but the spin-gap can, since for its determination only ground 
state energies of different spin sectors are needed. If inversion symmetry is present, one does 
not need to use complex phases and the FNMC method for a real wave function can be 
applied. 

In the present case, therefore, the FNMC appears to work well, in that the constraint 
imposed by the fixed node condition do not appear to have a dramatic effect on the energies 
and correlation functions over a reasonably large range of values of b. Unfortunately, as long 
as we lack more fundamental insight in the sign structure of the fermion wavefunction, it 
is difficult to say whether this is just one lucky example, or a robust property. Of course, 
in the present case our trial wave function has properly built in the tendency of the c and 
/-spins to form singlets, and our fixed-node estimates are not good in the extreme limits 
b^O (no local singlet correlations) and b— >1 (tightly bound singlets). 

For the smaller lattice sizes we have considered here, the sign problem does not appear 
to be so severeS, but for the larger lattice sizes, like those needed to study the spin triplet 
excitation or domain walls in de two-dimensional Hubbard mo dell, the advantages of the 
FNMC are more prominent. 

Our results also throw new light on our own earlier results for domain walls in the two- 
dimensional Hubbard mo del0. In these simulations, we compared ground state energies 
starting from a homogeneous trial state and from a domain wall trial state. Although the 
lowest energy state was found when applying a FNMC projection to a domain wall trial state, 
the ground state energy obtained after FNMC projection of a homogeneous trial state was 
found to be relatively close. Since it is conceivable that, e.g., domain wall type correlations 
do build up during the projection of a homogeneous trial state, a study of the correlation 
functions is needed before clear conclusions can be drawn. 
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